home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-02
/
tsptp.zip
/
FBENCH.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1993-04-09
|
13KB
|
379 lines
(******************************************************************************)
(* FBENCH.PAS *)
(* *)
(* This is a complete optical design raytracing algorithm, stripped of its *)
(* user interface and translated into "standard" Pascal. It not only *)
(* determines execution speed on an extremely floating point (including trig *)
(* function) intensive real-world application, it checks accuracy on an *)
(* algorithm that is exquisitely sensitive to errors. The performance of *)
(* this program is typically far more sensitive to changes in the efficiency *)
(* of the trigonometric library routines than the average floating point *)
(* program. This program has been demonstrated to be invariant under changes*)
(* in floating point format, as long as the format is a recognised double *)
(* precision format. If you encounter errors, please remember that they are *)
(* just as likely to be in the floating point editing library or the *)
(* trigonometric libraries as in the low level operator code. *)
(* *)
(******************************************************************************)
PROGRAM FBench(Output);
(******************************************************************************)
(* TIMING *)
(******************************************************************************)
(*$IFNDEF TopSpeed *)
(*%F TRUE *** Compile for Turbo Pascal ***)
USES TPBench;
(*%E*)
(*$ELSE *** Compile for TopSpeed Pascal ***)
IMPORT TSBench *;
(*$ENDIF *)
(******************************************************************************)
CONST
MaxSurfaces = 10;
VAR
CurrentSurfaces : BmInt;
Paraxial : BOOLEAN;
ClearAperture,
AberLSpher,
AberOsc,
AberLChrom,
MaxLSpher,
MaxOsc,
MaxLChrom,
RadiusOfCurvature,
ObjDist,
RayHt,
AxisSlopeAngle,
FromIndex,
ToIndex : BmReal;
SpectralLine : ARRAY [1..8] OF BmReal;
s : ARRAY [1..MaxSurfaces, 0..5] OF BmReal;
od_sa : ARRAY [FALSE..TRUE, 0..1] OF BmReal;
od_fline, od_cline: BmReal;
(** ASIN **********************************************************************
*
* Synopsis:
*
* PROCEDURE asin(x: BmReal): BmReal;
*
* Description:
* Returns the arc sine of x.
*
* NOTE:
* If x is out of rAnge 0.0 is returned but no further action is taken.
*)
FUNCTION asin(x: BmReal): BmReal;
BEGIN
IF (ABS(x) > 1.0) THEN
asin := 0.0
ELSE
asin := arctan(x / sqrt(1.0-x*x));
END;
(** COT ***********************************************************************
*
* Synopsis:
*
* PROCEDURE cot(x: BmReal): BmReal;
*
* Description:
* Resturns the cotangent of x.
* NOTE:
* If x is out of rAnge 0.0 is returned but no further action is taken.
*)
FUNCTION cot(x: BmReal): BmReal;
VAR
s: BmReal;
BEGIN
s := sin(x);
IF s = 0.0 THEN
cot := 0.0
ELSE
cot := cos(x) / s;
END;
(** TRANSITSURFACE ************************************************************
*
* Synopsis:
*
* PROCEDURE TransitSurface();
*
* Description:
* Calculate passage through a surface. If the variable Paraxial is TRUE,
* the trace through the surface will be done using Paraxial approximations.
* Otherwise, the normal trigonometric trace will be done.
*
* This routine takes the following inputs:
*
* RadiusOfCurvature Radius of curvature of surface being crossed. If 0,
* surface is plane.
*
* ObjDist Distance of object focus from lens vertex. If 0,
* incoming rays are parallel and the following must be
* specified:
*
* RayHt Height of ray from axis. Only relevant if
* ObjDist = 0.
*
* AxisSlopeAngle Angle incoming ray makes with axis at intercept.
*
* FromIndex Refractive index of medium being left.
*
* ToIndex Refractive index of medium being entered.
*
* The outputs are the following variables:
*
* ObjDist Distance from vertex to object focus after refraction.
*
* AxisSlopeAngle Angle incoming ray makes with axis at intercept
* after refraction.
*)
PROCEDURE TransitSurface;
VAR
iAng, (* Incidence angle *)
rAng, (* Refraction angle *)
iAngSin, (* Incidence angle sin *)
rAngSin, (* Refraction angle sin *)
oldAxisSlopeAngle,
sagitta: BmReal;
LABEL 999;
BEGIN
IF Paraxial THEN
BEGIN
IF RadiusOfCurvature <> 0.0 THEN
BEGIN
IF ObjDist = 0.0 THEN
BEGIN
AxisSlopeAngle := 0.0;
iAngSin := RayHt / RadiusOfCurvature;
END ELSE
iAngSin := ((ObjDist - RadiusOfCurvature) / RadiusOfCurvature) * AxisSlopeAngle;
rAngSin := (FromIndex / ToIndex) * iAngSin;
oldAxisSlopeAngle := AxisSlopeAngle;
AxisSlopeAngle := AxisSlopeAngle + iAngSin - rAngSin;
IF ObjDist <> 0.0 THEN
RayHt := ObjDist * oldAxisSlopeAngle;
ObjDist := RayHt / AxisSlopeAngle;
GOTO 999;
END;
ObjDist := ObjDist * (ToIndex / FromIndex);
AxisSlopeAngle := AxisSlopeAngle * (FromIndex / ToIndex);
GOTO 999;
END;
IF RadiusOfCurvature <> 0.0 THEN
BEGIN
IF ObjDist = 0.0 THEN
BEGIN
AxisSlopeAngle := 0.0;
iAngSin := RayHt / RadiusOfCurvature
END ELSE
iAngSin := ((ObjDist - RadiusOfCurvature) / RadiusOfCurvature) * sin(AxisSlopeAngle);
iAng := asin(iAngSin);
rAngSin := (FromIndex / ToIndex) * iAngSin;
oldAxisSlopeAngle := AxisSlopeAngle;
AxisSlopeAngle := AxisSlopeAngle + iAng - asin(rAngSin);
sagitta := sin((oldAxisSlopeAngle + iAng) / 2.0);
sagitta := 2.0 * RadiusOfCurvature * sagitta * sagitta;
ObjDist := ((RadiusOfCurvature * sin(oldAxisSlopeAngle + iAng)) *
cot(AxisSlopeAngle)) + sagitta;
GOTO 999;
END;
rAng := -asin((FromIndex / ToIndex) * sin(AxisSlopeAngle));
ObjDist := ObjDist * ((ToIndex * cos(-rAng)) / (FromIndex *
cos(AxisSlopeAngle)));
AxisSlopeAngle := -rAng;
999:
END;
(** TRACELINE *****************************************************************
*
* Synopsis:
*
* PROCEDURE TraceLine(line: BmInt; ray_h: BmReal);
*
* Description:
* Perform ray trace in specific spectral line.
*)
PROCEDURE TraceLine(line: BmInt; ray_h: BmReal);
VAR i: BmInt;
BEGIN
ObjDist := 0.0;
RayHt := ray_h;
FromIndex := 1.0;
FOR i := 1 TO CurrentSurfaces DO
BEGIN
RadiusOfCurvature := s[i][1];
ToIndex := s[i][2];
IF ToIndex > 1.0 THEN
ToIndex := ToIndex + ((SpectralLine[4] - SpectralLine[line]) /
(SpectralLine[3] - SpectralLine[6])) * ((s[i][2] - 1.0) /
s[i][3]);
TransitSurface;
FromIndex := ToIndex;
IF i < CurrentSurfaces THEN
ObjDist := ObjDist - s[i][4];
END
END;
BEGIN
WriteLn('FBench Benchmark.');
Paraxial := FALSE;
SpectralLine[1] := 7621.0; (* A *)
SpectralLine[2] := 6869.955; (* B *)
SpectralLine[3] := 6562.816; (* C *)
SpectralLine[4] := 5895.944; (* D *)
SpectralLine[5] := 5269.557; (* E *)
SpectralLine[6] := 4861.344; (* F *)
SpectralLine[7] := 4340.477; (* G *)
SpectralLine[8] := 3968.494; (* H *)
(*** Load test case into working array. The test case used in this ***)
(*** program is the design for a 4 inch achromatic telescope objective ***)
(*** used as the example in Wyld's classic work on ray tracing by hand, ***)
(*** given in 'Amateur Telescope Making', Volume 3. ***)
ClearAperture := 4.0;
CurrentSurfaces := 4;
s[1][1] := 27.05;
s[1][2] := 1.5137;
s[1][3] := 63.6;
s[1][4] := 0.52;
s[2][1] := -16.68;
s[2][2] := 1.0;
s[2][3] := 0.0;
s[2][4] := 0.138;
s[3][1] := -16.68;
s[3][2] := 1.6164;
s[3][3] := 36.7;
s[3][4] := 0.38;
s[4][1] := -78.1;
s[4][2] := 1.0;
s[4][3] := 0.0;
s[4][4] := 0.0;
(******************************************************************************)
(* Compute the looping overhead. The Dummy procedure must have some side- *)
(* effect so that it is not optimised out of existence. *)
(******************************************************************************)
StartTimer; (* Start the clock. *)
REPEAT
Dummy;
UNTIL NullTimesUp;
(******************************************************************************)
(* Now run the benchmark. Note that the Dummy procedure is also called so *)
(* that we can eliminate its overhead from the looping overhead. *)
(******************************************************************************)
StartTimer; (* Start the clock. *)
REPEAT
(*** Do main trace in D light. ***)
Paraxial := TRUE;
TraceLine(4, ClearAperture / 2.0);
od_sa[Paraxial][0] := ObjDist;
od_sa[Paraxial][1] := AxisSlopeAngle;
Paraxial := FALSE;
TraceLine(4, ClearAperture / 2.0);
od_sa[Paraxial][0] := ObjDist;
od_sa[Paraxial][1] := AxisSlopeAngle;
(*** Trace marginal ray in C. ***)
TraceLine(3, ClearAperture / 2.0);
od_cline := ObjDist;
(*** Trace marginal ray in F. ***)
TraceLine(6, ClearAperture / 2.0);
od_fline := ObjDist;
AberLSpher := od_sa[TRUE][0] - od_sa[FALSE][0];
AberOsc := 1.0 - (od_sa[TRUE][0] * od_sa[TRUE][1]) / (sin(od_sa[FALSE][1]) * od_sa[FALSE][0]);
AberLChrom := od_fline - od_cline;
MaxLSpher := sin(od_sa[FALSE][1]);
(*** D light. ***)
MaxLSpher := 0.0000926 / (MaxLSpher * MaxLSpher);
MaxOsc := 0.0025;
MaxLChrom := MaxLSpher;
Dummy
UNTIL BenchTimesUp;
(******************************************************************************)
ReportTimes;
(*** Reference results. These are derived from a run on Microsoft Quick ***)
(*** BASIC on the IBM PC/AT. ***)
WriteLn;
WriteLn('Reference results:');
WriteLn('Marginal ray: 47.09479120920 0.04178472683');
WriteLn('Paraxial ray: 47.08372160249 0.04177864821');
WriteLn('Longitudinal spherical aberration: -0.01106960671');
WriteLn(' (Maximum permissible): 0.05306749907');
WriteLn('Offense against sine condition (coma): 0.00008954761');
WriteLn(' (Maximum permissible): 0.00250000000');
WriteLn('Axial chromatic aberration: 0.00448229032');
WriteLn(' (Maximum permissible): 0.05306749907');
WriteLn;
(*** Calculated results. Any significant difference from the reference ***)
(*** results should be taken very seriously. This test has been ***)
(*** been demonstrated to be invariant for double precision floating ***)
(*** point formats. ***)
WriteLn('Calculated results:');
WriteLn('Marginal ray: ', od_sa[FALSE][0]:14, od_sa[FALSE][1]:14);
WriteLn('Paraxial ray: ', od_sa[TRUE][0]: 14, od_sa[TRUE][1]: 14);
WriteLn('Longitudinal spherical aberration: ', AberLSpher:14);
WriteLn(' (Maximum permissible): ', MaxLSpher:14);
WriteLn('Offense against sine condition (coma):', AberOsc:14);
WriteLn(' (Maximum permissible): ', MaxOsc:14);
WriteLn('Axial chromatic aberration: ', AberLChrom:14);
WriteLn(' (Maximum permissible): ', MaxLChrom:14);
END.